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Abstract: We analyse the statistical physics of a two dimensional lattice based gas 
with long range interactions. The particles interact in a way analogous to Queens 
on a chess board. The long range nature of the interaction gives the mathematics 
of the problem a simple geometric structure which simplifies both the analytic and 
numerical study of the system. We present some analytic calculations for the statics of 
the problem and also we perform Monte Carlo simulations which exhibit a dynamical 
transition between a high temperature liquid regime and a low temperature glassy 
regime exhibiting aging in the two time correlation functions. 
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1 Introduction 



It is believed that one of the essential ingredients for the formation of a glassy phase is 
the presence of frustration along with a large number of metastable states. The glass 
transition is a dynamical one, if one considers ordinary glass one sees the glass phase 
because on rapid cooling the system does not manage to find the crystalline structure 
of silicon dioxide. Instead the system becomes trapped on experimental time scales 
in a state of higher energy than the crystalline state. This phase is characterised by 
very slow dynamics and often exhibits aging phenomena in many observable quantities 
such as the correlation and response functions. A good introduction to the problem of 
the glass transition may be found in [|T|. Much of the understanding of glassy or out 
of equilibrium dynamics has been made via analogy with the dynamics of spin glasses 
and other systems with quenched disorder (for a review see [^). Ideas from the theory 
of spin glasses have been tentatively put forward to understand the glass transition ^ . 
In this paper we put forward a new model which we hope will be of use as a test of 
these ideas on systems without quenched disorder. We will emphasise its attractiveness 
as a model from both analytic and simulational points of view. 

A well known puzzle to chess players is the Eight Queens Problem. On a standard 
chess board it may be summarised as the problem of arranging eight queens such that 
no queen is in a position to take another queen at the next move. The interested reader 
will find that it is a non trivial task to find such an arrangement. In this paper we 
shall analyse a gas model based on the Eight Queens Problem. 

The attraction of this model is that the long range nature of the interactions facili- 
tates a rather straight forward and efficient implementation Monte Carlo dynamics for 
computer simulations. In addition we shall see that, whilst we cannot solve the statics 
exactly, various approximation schemes become analytically much simpler. 
The space of the model is an L x L lattice with periodic boundary conditions imposed 
to facilitate easier analytic and numerical analysis, the original problem of course does 
not have periodic boundary conditions. On this lattice are particles, the energy of 
a given configuration is given by the Hamiltonian 

N 

~ 2^ Xi,Xj "T ^yi,yj ~^ "xi-yi,Xj-yj ' "xi+yi,Xj+yjy \ ) 

where the position of the ith particle is denoted by the pair {xi,yi). The superscript 
L on the 6 indicates that they are the standard Kronecker delta functions but with 
arithmetic modulo L. One can see that for a periodic chess board the solution of 
the Eight Queens Problem is given the zero energy configurations of the Hamiltonian 
system above, the first two terms in the Hamiltonian represent the row and column 
constraints and the second two represent the constraint on the left to right and right 
to left diagonals respectively. 

In certain cases it is rather easy to find the zero energy states of the system. Let us 
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consider the case where L is prime. As an ansatz we shall take the configuration 

Xi = i and i/i = pi, (2) 

for each i between 1 and and where p is another a number different from L and 1. 
Clearly the column constraint is satisfied automatically. To violate the row constraint 
we would have to have 

p{i -j) = L (3) 

where = will always be taken to indicate equality modulo L, the fact that L is prime 
means that this may only be satisfied when i = j. Violation of the first diagonal 
constraint would mean that 

(p-l)(z-j) = L (4) 

This would imply that i = j or p — 1 = L. Violation of the second diagonal constraint 
would imply 

{p+m-j) = L, (5) 

hence either i = j or p+l = L. Therefore when one can find p such that p + 1 and 
p — 1 are not equal to zero modulo L then the above construction does indeed give a 
zero energy state. This analysis therefore demonstrates the existence of zero energy 
states for all L prime strictly greater than 3. 

This analysis is useful for our study of the dynamics; the ground states above are our 
analogy of the crystalline state of silicon dioxide when comparison is made with real 
glass. 



2 Mean Field Analysis 

One may carry out a geometrical mean field Boltzmann type of analysis to determine 
the energy per particle of the system. Consider the configuration about one particle 
at the point O which is fixed. Ignoring the effects of correlations the particle at O 
has two types of sites surrounding it, those with which it interacts directly, i.e. those 
which are on the same row, column or diagonal which we shall call type A, and those 
with which it does not interact directly which we shall call type B. Denote by Nx^y 
the number of sites of type Y that interact with a given particle of on a site of type X 
{X, Y e {A, B,0}). Simple counting yields 

No-^A — 1, by definition 

Na-^a — N + 7, interaction with A" — 2 on the line interacting with O 

and intersection with 3 A s on the 3 other lines 

Nb^a = 3A^ - 12, using Nq-^a + Na^a + Nb-.a = 4Ar - 4 

Nq-^b — 0, by definition 

Na-^b — 12, 4 lines with 3 intersections each 

Nb-.b = m - 16, using No-.a + Na-.a + Nb-.a^4:N- A (6) 



2 



The Boltzmann equations for the average number of particles in sites of type type A 
and B are 

(A) = Aexp (-(TV + 7)P{A) - (37V - 12)P{B) - /3) (7) 

and 

(B) = X exp i-12l3{A) - {AN - 16)P{B)) . (8) 

Here A is a Lagrange multipher factor enforcing the over all particle number to be N 
i.e so that 

1 + A(N-1){A) + {N^-AN + 3){B)=N, (9) 

which simplifies to 

4:{A) + {N-S){B)^1. (10) 

In the limit TV — > oo one finds 

{B) = 1/N and {A) = a/N, (11) 
where a is finite. Consequently 

(B) = A(l-^(l-e-/^))-x(l-l(l-e-0)^^-^ 

Aexp (-4(1 -e" 
_ 1 
N' 

Hence 



N 

Taking similarly the large N limit for the equation for (A) yields 



a — exp 



(12) 



A=lexp(-4(l-e-^)). (13) 



:i-e-^)(l-a)-/3), (14) 

and the energy E per particle is given by 

E{P) = 2a{P). (15) 

For infinite temperature one finds that E{0) = 2 and for low temperatures E{f3) ~ 
2e^~^. Hence the mean field calculation does not contradict the fact that there are 
zero energy ground states for N prime. 
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3 The Hypernetted Chain Approximation 

Perhaps one of the most successful approximations in gas and hquid theory is the hy- 
pernetted chain (HNC) approximation , which is an integral equation which resums 
a certain class of diagrams in the virial expansion. Normally this integral equation 
must be resolved numerically and its resolution is rather difficult. However we shall 
see here that for our model the resulting equation is drastically simplified in the large 
N limit. 

The matrix form of the HNC equation is 

log(l + hij) = -pVij +pY, hikihkj - log(l + hkj) - (3Vkj) (16) 

k 

where ^ 

= ^ ( ( APi ) - P^ii - ) • (17) 

Here pi is the density at the site i and p = (p) in homogeneous systems. In the present 
model p = 1/N. One may write the interaction matrix as 

= 35,, + Wij (18) 

The important simplification comes from noting that (using the counting from the 
previous section) 

W"^ = {"iN -9)1 + {N -Q)W + 12U, (19) 
where / is the identity matrix and U is the matrix with each element equal to 1. Also 

U"^ = N'^U and UW = WU = {AN - 3)U. (20) 

Hence the elements W, I, U, form a closed algebra. In general we may therefore repre- 
sent an element in this algebra by a triple 

A = (ai, a2, 03) = asU + (02 - a3)W + (ai - 02)/ (21) 

where 

Aij = ai when i = j 

^ij = 0'2 when i ^ j and Vij = 1 

Aij = 03 when Vij = 0. (22) 

In this notation the product of two elements A = (01,02,03) and B = (6i,&2;^3) is 
given by C = AB = (ci, C2, C3) where 

ci = 03&3 A^' + («i - 02) ih - 62) + (a2 - 03) (62 - h) {AN - 3) 

+ (03(^2 - &3) + h{a2 - a3)){AN - 3) + (03(61 - 62) + h{ai - 02)) 

+ {{ai - a2){b2 - h) + {h - b2){a2 - a^)) (23) 
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C2 = ashN^ + {a^ - a3){b2 - b3){N + 6) 

+ (03(^2 - &3) + h{a2 - a3))(4iV - 3) + (03(61 - 62) + h{ai - as)) 

+ {{ai - a2){b2 - ba) + {bi - b2){a2 - a^)) (24) 



and 



C3 = a363A^' + 12(a2-a3)(62-63) 

+ ia3{b2-bs) + b3ia2-as)){4:N -3) + {as{bi-b2) + b3{ai-a2)) (25) 

Using this algebra and taking the hmit — 00 at the end of the calculation one finds 
that in this representation the HNC reduces to 

i + ft2 — ^3 

log{l + h2) = -(3 + —— — , 27 

1 + A12 - AI3 

and 

log(l + h) = h. (28) 
The final equation is trivial giving = and hence 

h, = exp{-4p + - 1 (29) 

1 + /i2 

with 

log(l + M = -/3 + -4V- (30) 

1 + /l2 

For /3 small one finds /12 P and for /3 large /la l+l/(/?+2-log(/?)+0(log(/?)//?)). 

The energy per particle is then given by 

N 

^ = T E(PoP.)V^)i = 2(/i2 + 1). (31) 
^ i=i 

Hence, at high temperatures, E ~ 2(1-/5) as in the mean field calculation, but 
E ^2/{P + 2- \og{/3) + 0{\og{l3)/l3)) at low temperatures. 



4 Monte Carlo Simulation 

The dynamics for the Hamiltonian (|ip were simulated using a Monte Carlo method 
using the Metropolis algorithm a with sequential update on the particles. Two variants 
of the dynamics have been analysed, firstly a non local dynamics where a particle may 
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move to any site on the board and secondly a local random walk dynamics where a 
particle may move to any of its eight nearest neighbours on the lattice. One of the 
advantages of the current model is that the energy may be expressed as a function of 
the occupation numbers of each row {Nr{k)), column {Nc{k)), left to right diagonal 
{Nd-{k)) and right to left diagonal {Nd+{k)). This speeds up greatly the calculation 
of the energy change at each move. In this notation it is easy to see that 

1 ^ 

^ = ^ T^i^rik)' + N,{kf + Nd-{kf + iV,+ (A;)) - 2 (32) 

k=l 

This leads to a reduction by a factor of of the time needed to compute the energy 
change with respect to a calculation using the particle positions. The form of the 
energy given by equation (^) demonstrates that if one were to neglect the diagonal 
interactions one would arrive at a model similar to two independent Backgammon 
models but at negative temperature. The Hamiltonian is not strictly the same 
but the tendency is to put all the particles in different boxes rather than to put them 
into single box which is the ground state of the Backgammon model. Here the free 
energy barrier clearly have an energetic component as opposed to the Backgammon 
model where their origin is entropic. We shall also see later that, in common with the 
Backgammon model, the dynamical transition temperature appears to be at T = 0. 
The system was started from a random initial configuration, thus simulating a rapid 
quench from high temperature to the simulation temperature. In both the cases of non 
local and local dynamics the energy as a function of time was measured. 
In the case of non local dynamics the particle-particle correlation function treating the 
particles as distinguishable. It is defined by 

i 

note that C is normalised to be one when t = t' (note in keeping with convention t' 
will always denote the earlier time). In the case of local dynamics it is interesting to 
measure the correlation function which gives us information about the effective particle 
diffusivities, this correlation function is normally denoted by B[t,t') in the literature 
and is given by 

m') = ^Ei^^it)-^^it')r. (34) 

i 

We shall discuss the results of the two dynamic types separately. 
Non Local Dynamics 

With non local dynamics one expects, at an intuitive level, that the dynamics should be 
faster and less prone to be glassy with respect to local dynamics. If one imagines a self 
consistent picture of a non ordered phase, then a single particle moves in a potential 
generated by the other particles. With non local moves trapping mechanisms at the 
spatial level are less important. 
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The simulations we carried out were for systems of 50000 particles and times of up to 
60000 Monte Carlo sweeps through the entire system. In addition we verified that for 
systems of size 48611 (which is prime) the results were not significantly altered. 
In figure (1) is shown the results of the energy per particle measured by the simulation 
for t = 30000. Up to values of /3 = 5 we found that hj t = 30000 the energy density 
had reached a well defined asymptotic value. We found for values of /3 > 6 that at 
t = 30000 the system had not reached an asymptotic values of the energy density and 
that the value of the energy density continued to decay though very slowly. The values 
of E shown on figure (1) for /3 > 6 are the values measured at t = 30000 and are 
not equilibrium values. One sees from figure (1) that between (3 = 5 and /3 = 6 the 
measured energy flattens off rather suddenly in a fashion reminiscent of glassy systems. 
For vales of /? < 5 and a minimum waiting time of 10000 we found equilibrium behaviour 
for the correlation function. In this region 

C{t,t') = C{t-t') ^ Aexp{^-^). (35) 

To 

One may regard Tq clS cl characteristic time scale, if this time scale is greater than 
the observation time then one should expect out of equilibrium behaviour. Shown 
in figure (2) is a plot of log(ro(/3)) as measured by fitting the simulation data with 
the formula (|35D above. As one can see the plot is linear and the best fit is To{(3) = 
0.317 ■ exp(l. 796/3), when /3 is such that ro(/3) is much less than the age of the system 
one expects that equilibrium should be achieved. For an age of the system of order 
10000 one finds that the value of the temperature where the age and ro(/3) have the 
same age is /3 ~ 5.7 which is entirely consistent with the plateau observed in the energy 
and with the breakdown of time translational invariance in the correlation function. 
For values of /3 > 6 we found that the function C{t, t') was no longer time translational 
invariant, neither did it have an exponential form. For values of /5 > 7 we found that to 
a reasonable degree of accuracy the system exhibits perfect aging, i.e. C{t, t') = f{t/t'), 
see for example the curves in figure (3) for (3 = 10. At T = one sees that the 
system also ages and the five curves shown collapse onto the same master curve when 
plotted against 7^, this is shown in figure (4). This is presumably because at T = 
the characteristic time tq becomes infinite and there is no interrupted aging p. The 
scaling t/t' is been observed in a variety of mean field models and also in the 
phenomenological trap model P]. The improvement of the perfect aging scaling as one 
approaches zero temperature was also observed by Ritort in Monte Carlo simulations 
of the Backgammon model 0] . 
Local Dynamics 

Here we consider the more physical situation of local random walk dynamics. 
Shown in figure (5) is the energy measured for the local dynamics for a system of size 
50000 particles up till a time of 2000 Monte Carlo sweeps. Up till /3 = 4 the energy 
reaches an asymptotic value, where as for (3 > 4 the energy continues to decrease. The 
values shown in these cases is the minimum value of the energy attained. As one sees 
clearly in the figure, this measured value is not monotonic and actually increases after 
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(3 — 6.5. For longer simulation times this minimum in the measured dynamic energy 
shifts towards larger vahics of 

For values of the waiting time t' = 10000, 20000, 30000, 40000 up till /3 = 4 one sees 
that 

B{t,t') = K{P){t-t') (36) 

and hence we are in an equilibrium phase for these temperatures and at this scale of 
age. The values of \og{K,{/3)) are shown in figure (6) and one sees that the curve is 
linear and one finds that k,{P) — 0.34 • exp(— 1.08/3). 

In the nonequilibrium regime we find that one may fit the curves as an anomalous 
diffusion but with a prefactor that depends on the waiting time. The behaviour of the 
particles is sub diffusive and we have fitted it with the form 

B{t,t')^c{P)^-^^, (37) 

where 7 = 2(1 — a). This simple linear relation between the exponents for {t — t') and t' 
works rather well but we have no theoretical reason to expect this result. An example 
of the rescaled curves for /3 = 6 is shown in figure (6). The values of a as a function of 
/3 are shown in figure (7). To go to lower values of the temperature than those shown 
is rather difficult as the noise in the measured values of B{t,t') becomes rather large, 
presumably because one is approaching the zero temperature transition. 



5 Conclusions 

We have presented a preliminary study of a two dimensional particle system which 
exhibits a dynamical transition between a gas/liquid and a glassy phase. As expected 
from the dynamical nature of the transition the temperature at which one sees this 
transition is dependent on both the time scale of the observations and on the nature of 
the dynamics. The long range nature of the interactions leads to an compact application 
of various standard approximation schemes for the statics - particularly the hypernetted 
chain approximation which works reasonably well in the high temperature phase. The 
model is also attractive from the simulation point of view as the form of the interactions 
leads to an efficient implication of the Monte Carlo updating procedure. 
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List of Figure Captions 

Fig 1. Nonlocal dynamics: Dynamic Energies for L = 50000 at t = 30000. Calculated 
static energies using the HNC and the mean field approximation are also show. 
Fig 2. Nonlocal dynamics: logro(/3) L = 50000 (with linear fit shown). 

Fig 3. Nonlocal dynamics: C{t, t^) plotted against t/t^ at /5 = 10 for t^ = 5000, 10000, 15000, 20000, 2 
(L = 50000). 

Fig 4. Nonlocal dynamics: C{t,ty,) plotted against t/t^ atT = Ofort^ = 20000, 40000, 60000, 80000, 1 
(L = 50000). 

Fig 5. Local dynamics: Dynamic Energies for L = 50000 at t = 2000. 
Fig 6. Local dynamics: —\og{K{j3) in the equilibrium regime (with linear fit shown). 
Fig 7. Local dynamics: Rescaled B{t + t^,t^) for t^ = 10000, 20000, 30000, 40000 
at (3 = 6 with a taken to be 0.862. 

Fig 8. Local dynamics: a(/5) in the nonequilibrium regime. 
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